Beyond the local approximation to exchange and correlation: the role of the Laplacian 

of the density in the energy density of Si. 
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We model the exchange-correlation (XC) energy density of the Si crystal and atom as calculated 
by variational Monte Carlo (VMC) methods with a gradient analysis beyond the local density 
approximation (LDA). We find the Laplacian of the density to be an excellent predictor of the 
discrepancy between VMC and LDA energy densities in each system. A simple Laplacian-based 
correction to the LDA energy density is developed by means of a least square fit to the VMC XC 
energy density for the crystal, which fits the homogeneous electron gas and Si atom without further 
effort. 

PACS numbers: 71.15.Mb, 31.15.Ew, 71.10.Ca 



The crucial ingredient of density functional theory P, 
2] (DFT) is the the exchange-correlation (XC) energy 
which incorporates the effects of many-body correlations 
on the ground-state energy of an electronic system into 
its expression as a functional of the ground-state den- 
sity. The success and widespread application of DFT 
in solid-state physics and quantum chemistry has been 
due to the remarkable accuracy of simple and efficient 
local and "semilocal" models for this quantity, including 
the local density approximation (LDA) N] ; generalized 
gradient approximations (GGA's) 0, 0, U j an d various 
extensions of the GGA [a, 0, @- These methods form 
a hierarchy of approximations in which this intrinsically 
nonlocal and as yet poorly understood functional of the 
density is mapped to a succession of increasingly complex 
local functions of the density, its gradient and related 
quantities. However, no systematic method for develop- 
ing such corrections is known to exist, and the accuracy 
of current methods is not yet consistently at the level 
(roughly a milli-Rydberg) needed to characterize chem- 
ical reactions and other applications highly sensitive to 
the total energy. 

A fruitful source of intuition and of mathematical con- 
straints in the development of DFT's has been the analy- 
sis of the XC energy in terms of the XC hole, the change 
in density from the mean that occurs about an electron's 
position due to exchange and Coulomb correlations Q. 
It provides a natural interpretation for the XC energy 
density and thus has aided in the construction of several 
DFT models 0, @, @ Despite the usefulness of the XC 
hole in DFT development, there have been few calcula- 
tions of it for realistic systems. Recently, however, accu- 
rate variational Monte Carlo (VMC) calculations of the 
XC hole and the associated e nerg y density have been per- 
formed for the Si crystal 0, Hlf and atom [l^ within a 
pseudopotential approximation. These calculations have 
provided a wealth of data for analysis jlflj , but a compre- 
hensive understanding of their implication for DFT has 
to date been lacking. 



We present in this paper an analysis of the XC energy 
density associated with the XC hole in the Si crystal and 
atom in terms of a gradient analysis of the density. We 
find that the deviation of the XC energy density from the 
LDA model is markedly correlated with the local Lapla- 
cian of the density, a quantity that has been mostly ne- 
glected in developing DFT's, with the local gradient play- 
ing little or no role. We construct a minimal Laplacian- 
based model to quantify this relation with parameters fit 
to the crystal data. The resulting fit captures most of 
the discrepancy between the VMC and LDA energy den- 
sities, and fits both the homogeneous electron gas (HEG) 
and Si atom cases with no further effort. 

A strong correlation between the Laplacian of the den- 
sity and the XC energy density has previously been re- 
ported ^3 f° r a model strongly inhomogeneous electron 
gas. However, the current work is the first time that such 
a picture has been found in the context of the complex- 
ities (covalent bonding, atomic orbitals, diamond struc- 
ture) inherent in a real material, one that is paradigmatic 
for all covalently bonded systems. These results suggest 
the existence of a simple yet universal correlation be- 
tween the XC hole and the local density Laplacian that 
should be a help in guiding future DFT models. 

In DFT, the XC energy E xc is usually written as an 
integral of a locally defined XC energy density, e xc : 



Er 



d 3 r e xc {r; [n]); 



(1) 



where e xc is itself an unknown functional of the density 
n. The simplest ansatz for e xc is that of the LDA in 
which the true nonlocal functional at a given point in 
space is replaced by that of the homogeneous electron gas 



(HEG) with the local value of the density: 



.LDA 



(r; [„]) 



e xc EG (r s (r)), where r s = (3/47m) 1/3 is the Wigner-Seitz 
radius. Corrections to the LDA are usually based on 
a gradient expansion in which the variation in the 
density near r, described by derivatives of n(r), is used 
to modify e xc (r). GGA's add a dependence on |Vn(r)|, 
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FIG. 1: Comparison of DFT and VMC XC energy densities on the (110) plane of the Si crystal, (a) Difference in the LDA 
XC energy density and that of VMC data [Tol lll|. (b) Difference between that of the "GGA ++ " model described in the text 
and VMC. Contours in increments of 2xl0 -3 a.u., with thicker contour that for zero energy difference. Bluer regions show 
negative difference and redder regions, positive. 
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FIG. 2: Gradient analysis of the density of crystalline Si. The density (a), its gradient squared (b), and Laplacian (c) on 
the (110) plane of the Si crystal. Atoms and bonds outlined in black. Shading varies from blue (lowest) to red (highest) and 
contours are in increments of 0.01 a.u. (a), 0.01 a.u. (b), and 0.05 a.u. (c). In (c) the zero contour is the thicker black line. 



and mctaGGA's |7|, on more complex local derivatives. 
The density Laplacian V 2 n occurs to the same order as 
|Vn| in the gradient expansion but is less often used 0. 

An intuitive picture of e xc is obtained from the XC 
hole n xc (r, r'), which measures the change in density at r' 
from the mean density n-(r'), given the observation of an 
electron at r. The XC energy density may be expressed 
in terms of the adiabatically integrated XC hole jig : 

e xc (r) = ^fdX fdh'^l. (2) 
1 J Jo l r ~ r 

(In this paper, we use hartree atomic units.) Here, 
n xc represents the XC hole evaluated for a system with 
Coulomb coupling Ae 2 and the same ground-state density 
n(r) as the true system. In this formalism, ea; C (r)/n(r) 
is the sum of the potential energy due to the interaction 
of an electron with its own hole and the kinetic energy 
cost to create the hole. 

Unfortunately, e xc is not uniquely definable - any func- 
tion that integrates to zero over the system volume could 
be added to e xc in Eq £Q to generate a new "gauge" 
choice for the energy density, to which the energetically 
relevant quantity E xc would be invariant. This is implic- 
itly done in GGA's to convert any potential dependence 



of e xc upon V 2 n to an equivalent dependence upon Vn| 2 
alone 15]. On the other hand, the adiabatic method is a 
natural, easily interpreted choice for defining e xc ; more- 
over it is readily calculable in the VMC method from 
the expectation of the XC hole taken for several different 
values of A flfij . 

To visualize the task faced in describing e xc for a re- 
alistic system, we plot in Fig. ^a) the difference 5e xc 
between the e xc of the LDA and that of the VMC calcu- 
lation of Hood et al. 0, El for the Si crystal in a pseu- 
dopotential approximation. The LDA predicts too deep 
an energy in the region of the Si bond, and too shallow 
an energy at low density, most obviously in the pseudo- 
atom core, but also, amplified in effect since it includes a 
large percentage of the unit cell volume, in the interstitial 
regions of the crystal. The net contribution of positive 
and negative errors in e xc almost exactly cancel, so that 
the integrated E xc in the LDA is essentially the same as 
that of the VMC [ll| . The exact functional behavior of 
the energy density difference is quite complex. 

Figure [3 shows a gradient analysis of the Si-crystal va- 
lence electron density on the (110) plane. The gradient of 
the density squared |Vn| 2 , shown in Fig[2Jb), highlights 
the critical points of the density as blue regions where the 
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gradient is nearly zero. It is significantly nonzero around 
the edges of the bond between two Si atoms. The Lapla- 
cian V 2 n, Fig. |2fc), is negative in regions of strong elec- 
tron localization in the bond and positive in regions of 
electron depletion, such as the atom core and the inter- 
stitial regions. It has a characteristic "butterfly shape" 
in the bond center, caused by two regions of peak density 
located near the two Si atom valence shell peaks. 

Upon comparison of Figs.^a) and [2 what is immedi- 
ately evident is that the shape delineated by V 2 n charac- 
terizes the discrepancy between the VMC and LDA XC 
energies. It reliably predicts the sign of the correction 
needed on a point by point basis throughout the unit 
cell, identifies regions of maximum error (bond and atom 
core), and reproduces key topographic features such as 
the shape of the region of maximum energy error in the 
bond. In contrast |Vn| 2 seems to have little to do with 
the trends in energy density error. 

VMC calculations of e xc have recently been performed 
for the valence shell of the Si atom in a pseudopotential 
model [13 ■ These allow us to verify the trends demon- 
strated in the gradient analysis of the crystal in a system 
that lacks bonds, and has significantly different boundary 
conditions. Shown in Fig. Of a) are |Vn| 2 and V 2 n of the 
Si pseudo-atom electron density versus radial distance. 
The peak of the density, indicated by the vertical dotted 
line, marks the zero of Vn and the maximum negative 
value of V 2 n. The solid line in 3(b) shows the differ- 
ence in e xc between the local spin density (LSD) [Hi 
and VMC results. Ignoring short-wavelength statistical 
fluctuations that are a by-product of the Monte Carlo 
calculation, a dramatic correlation of Se xc with V n is 
seen, with the same qualitative trends as the crystal. 

These two examples (crystal and atom) demonstrate 
a qualitatively consistent dependence of e xc upon the 
Laplacian of the density that should be quantifiable - 
but not in the context of GGA's, which do not include 
a dependence on V 2 n. We consider an enhanced GGA 
model, a "GGA++", of the form 



e^++(r s , S 2 



I) = F xc {r s 



(3) 



where the correction to the LDA energy density is 
expressed by an enhancement factor F xc dependent 
upon the Wigncr radius r s and dimensionless variables 
I = r 2 (r)V 2 n(r)/n(r) and s = r s (r) |Vn(r)| /n(r). This 
GGA++ is fit to VMC data for the Si crystal by a least- 
squares procedure that minimizes the variance in the en- 
ergy density from the VMC value, integrated over the 
unit cell. The root-mean-square error of the energy den- 
sity, 5e rms , obtained in this way is 0.442 millihartrees for 
the LDA, and represents the average deviation from zero 
for the energy-difference plot shown in Fig. ^a). 

We have found that a form for F xc depending only 
upon the dimensionless Laplacian / provides the optimal 
fit to our data in the sense of returning the greatest de- 
gree of correction per fitting parameter. The form is 
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FIG. 3: Gradient analysis and GGA ++ fit of e xc for the 
Si pseudo-atom, (a) Gradient squared and Laplacian of the 
density as a function of radial distance from the atom core, 
(b) Difference between e xc of the LSD and GGA ++ models 
and that obtained from VMC p^| . 



given by 



F XC {1) 



1 



a + (31 



(4) 



with optimized fitting parameters a = —0.0007, (3 = 
0.0080, and 7 = 0.026. The fitting error 5e rms is thereby 
reduced 70% from its LDA value to 0.132 millihartrees. 
This form potentially satisfies several known properties 
of the universal e XCl particularly recovering the correct 
value in the HEG limit (s 2 = 1 = 0) for a = 0. It behaves 
properly under uniform scaling to infinite density [l8j but 
fails to include a dependence of F xc on r s due to correla- 
tion. The smallness of the optimized value of a indicates 
that the best fit for the Si crystal simultaneously satisfies 
the HEG limit. This supports the validity of our model 
as a description of a genuine physical phenomenon rather 
than a mathematical anomaly specific to Si. 

Shown in Fig. ^b) is the difference in energy density 
between our three-parameter GGA ++ fit and the VMC 
data of Hood et al. , on the same energy scale as the en- 
ergy difference between LDA and VMC in (a), showing 
point by point what 5e rms shows on average. The differ- 
ence in e xc has been greatly reduced everywhere through- 
out the unit cell, with the exception of the bond center 
and at the antibond point behind each bond. We have 
also tried a 5-parameter fit including terms of order s 2 
and forms with higher order corrections, with only mini- 
mal improvement of 5e rms . In every case tried the linear 
coefficient (3 for I remains at 0.008 to within 10%. 

The transferability of our model can be tested by ap- 
plying it to the Si atom data of Ref. 0. We have ap- 
plied the Laplacian-only F xc obtained from our fit to the 
crystal data without any further adjustments as a cor- 
rection to the LSD XC energy density for the Si atom. 
This is defined similarly to Eq. 0, by e ^SD-GGA++ = 



4 



F xc (s 2 ,l)e^ D (r s ,0, where C 



is the local spin 



polarization. The result for 5e xc using this model is 
shown in Fig. E^b) ; the overall error 5e rms is reduced by 
70% from its LSD value, achieving the same reduction of 
error as for the crystal. 

Our numerical results tieing e xc to V 2 n can be mo- 
tivated qualitatively by reconsidering a gradient expan- 
sion, this time for n xc . This would use as input the 
change in density within the length-scale of the XC hole 
about any position, as described by local derivatives of 
the density, to correct the errors inherent in the LDA as- 
sumption of a locally homogeneous environment. As the 
Coulomb interaction is dircctionally invariant, only the 
change in density averaged over angle should contribute 
to this correction. This is precisely what is measured 
by V 2 n, and is unobtainable from |Vn| 2 . Given an e xc 
derived from the adiabatically integrated XC hole, one 
could then expect the error in the LDA model of e xc to 
be dominated by the local value of V 2 n. 

The value of the Laplacian of the density in electronic 
structure has been noted in several other contexts. It has 
been used successfully as a diagnostic tool in character- 
izing the electronic structure of molecules 0] . Covalent 
bonds have been found to be distinguished by a negative 
Laplacian at the bond center, denoting the build-up of 
charge within the bond, and non-covalcnt ones by a pos- 
itive V 2 n; in addition the hour-glass pattern observed 
in the Si crystal bond is typical of other tetrahedrally 
bonded systems. Secondly, studies of the XC potential 
of atoms [13, El have pointed out that terms in V 2 n 
are necessary to model the potential in the nuclear cusp 
and asymptotic regions. Thus the relevance of this quan- 
tity to DFT extends beyond the pseudopotcntial mod- 
els studied here to all-electron calculations, and possibly 
from covalent to other types of chemical bonds. 

The XC potential v xc (r) = 8E XC /Sn(r), necessary for a 
self-consistent determination of the density is easily ob- 
tained within the plane-wave pseudopotential formalism 
of the DFT. Self-consistent calculations of density and 
structural properties of Si using our GGA ++ model show 
no significant deviation from the already reasonably good 
prediction of these quantities in the LDA. Full results will 
be discussed in a further paper. 

A caveat in regard to our results is that our model has 
been fit to data obtained by a variational method that 
underestimates the correlation energy. The true correla- 
tion energy for each system may be lower than that of 
the VMC by about 15%, and E xc lower by 1-2%. How- 
ever, within the VMC approximation, the main effect of 
adding correlation has been to increase the match be- 
tween the LDA error and V 2 n from that observed in the 
exchange-only case, shown in Fig. 6(a) of Ref. El The 
effect of the addition of the missing correlation energy 
might well be to reduce further the discrepancy between 
the actual e xc and a Laplacian fit. 

In summary, our fit of e xc in terms of a Laplacian based 



enhancement factor F XC {1) provides a simple model that 
has a surprisingly wide range of applicability: from the 
HEG to covalently bonded crystal to open shell atom. 
This points to the potential for a Laplacian-based F xc to 
make an excellent approximation to the true, universal 
one for a wide range of systems. To date, the devel- 
opment of GGA's and metaGGA's has emphasized the 
gradient of the density as the basic departure point for 
the post-LDA description of DFT. Our analysis indicates 
rather that it may be advantageous to start with V 2 n as 
the key factor in going beyond the LDA. 
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